home *** CD-ROM | disk | FTP | other *** search
/ Precision Software Appli…tions Silver Collection 4 / Precision Software Applications Silver Collection Volume 4 (1993).iso / stats / chadyn.exe / YLYAP.C < prev    next >
C/C++ Source or Header  |  1988-11-28  |  22KB  |  744 lines
  1.  
  2. /*******  YLYAP.C  for computing Lyapunov numbers and NEWTON's Method  *******/
  3. /*********************** (C) 1986,7,8 by JAMES A. YORKE **********************/
  4.  
  5.  
  6. #define EL if(printer >=1) erase_line()  
  7.  
  8. /***********     Routines in YLYAP      *************
  9.  
  10. lyapinit()  this should be done when dot = -preiter but isn't 
  11. lyapunov()  one step in computing lyap numbers corresponding to one 
  12.                 step of the differential equation;
  13.     this is called by iterate_map() 
  14. print_lyap(output)
  15. set_y()   y[j] is set for j = 5,...,20   -- for dif eq 223 
  16. set_vector()
  17. orthonormalize()  orthonormalize  the vectors 
  18. orthogonalize(i)  change vector[][i] into vector[][i] so that vector[][i] 
  19.     is perpendicular to previous ones. This new vector[][i] is 
  20.     unnormallized 
  21. double dot_p(kk,ii)  returns dot product
  22.  
  23.   Newton Routines 
  24. newtinit()
  25. newton(N)
  26. ************************************************************/
  27.  
  28.  
  29.  
  30. /* vec_dim is the number of coordinates  of the Lyapunov vectors 
  31.     so jj satisfies 0 <= jj < vec_dim    */
  32. /* num_lyap is the number of Lyapunov exponents to be computed. 
  33.     The default value is 0  and the max is vec_dim     
  34.     so i satisfies 0 <= i < num_lyap */
  35. /* lyapzero : y[lyapzero] is the zeroth coord of the zeroth Lyapunov 
  36.     vector */
  37. /* these three are all set by select_map() in YMAPS.C. 
  38.     In addition that routine sets init_process() to be something like 
  39.     jac_223() to set appropriate
  40.     constant jacobian terms in Ydes */
  41.  
  42. #include "yinclud.h"
  43. static double   vector[5][5],
  44.                 dot_prod[5];
  45. static double   product;    /* this is for computing products of norms of
  46.                    the vector[][] vectors   */
  47.  
  48. /* in 223, y[0] is time,y[1] and y[3] are x variables 
  49.     and y[2] and y[4] are x-primes */
  50.  
  51.  
  52. lyapinit() {            /* this should be done when dot = -preiter and
  53.                    when num_lyap is changed especially if it is
  54.                    done while the program is iterating */
  55.     int    i;
  56.     int    coord;
  57.  
  58.     lyaptime = 0;
  59.     if(step != -9999) {
  60.         lyapstep = step * its_per_plot;/* dif equation case */
  61.         dim = lyapzero + num_lyap * vec_dim;/* for rungekutta */
  62.     }
  63.     else
  64.         lyapstep = its_per_plot;/* map case */
  65.     for(i = 0; i < num_lyap; i++)
  66.         for(coord = 0; coord < vec_dim; coord++)
  67.             y[lyapzero + coord + vec_dim * i] = 1.+ coord + vec_dim * i;
  68.  
  69.     set_vector();
  70.     orthonormalize();
  71.     set_y();
  72.  
  73.     for(i = 0; i < vec_dim; i++)/*  initialize vector lyapsum  */
  74.         lyapsum[i] = 0;
  75.     if(level >= PROCESS && printer > 0) {
  76.         scr_rowcol(1, 0);
  77.         PRINT "\nLyapunov exponents are now(re-)initialized.");
  78.     }
  79. }
  80.  
  81.  
  82. lyapunov() {            /* This routine takes the Lyapunov vectors and
  83.                    orthonormalizes them; in the process, it
  84.                    computes growth and shrinkage rates;
  85.                    sometimes it prints out estimates of
  86.                    exponents; one step in computing lyap
  87.                    numbers corresponding to one step of the
  88.                    differential equation; it is called by
  89.                    processi() in Q.C */
  90.     FILE * output;
  91.     int     puts();
  92.     output = StOutPut;    /* this step of introducing output rather than
  93.                    using StOutPut directly seems stupid but
  94.                    print_lyap(stdout) doesn't seem to work with
  95.                    the big model DeSmet */
  96.     if(dot == 0) {
  97.         lyapinit();    /* this should really be done at the outset of
  98.                    computation */
  99.     }
  100.     else {
  101.         set_vector();
  102.         orthonormalize();
  103.         set_y();
  104.         if(printer > 2) {
  105.             scr_rowcol(0, 0);
  106.             print_lyap(output);
  107.             if((10 * (dot / 10)) == dot) {
  108.                 /* print warning evey tenth dot */
  109.                 erase_line();
  110.                 PRINT
  111.                     "Hit 2 to stop printing & speed computing lots\n");
  112.             }
  113.         }
  114.     }
  115. }
  116.  
  117. print_lyap(output)
  118. FILE * output;
  119. {
  120.     int     num = num_lyap;
  121.     int    i;
  122.  
  123.     if(num_lyap < storeNumLyap)
  124.         num = storeNumLyap;/* needed for straddle orbits */
  125.     if(lyaptime > 0) {
  126.         if(level >= PROCESS) {
  127.             fprintf(output, "\n");
  128.             EL;
  129.             fprintf(output, "\n");
  130.             EL;
  131.             fprintf(output, "dot =%ld; ", dot);
  132.                 /* no '\n' wanted */
  133.         }
  134.  
  135.         EL;
  136.         fprintf(output, "Lyap exps= ");
  137.  
  138.         for(i = 0; i < num; i++)
  139.             fprintf(output, "%lf ", lyapsum[i] / lyaptime);
  140.         lyap_dim(output);
  141.     }
  142.     else
  143.         if(lyaptime < 0)
  144.             fprintf(output,
  145.                     "No lyap exponent totals computed yet; lyaptime = 0             \r");
  146.         else
  147.             if(lyaptime == 0 && level >= PROCESS
  148.                     && output == StOutPut)
  149.                 /* for clearing out old print on screen */
  150.                 erase_line();
  151. }
  152.  
  153. lyap_dim(output)
  154. FILE * output;
  155. {
  156.     int     i;
  157.     double  SumOfExponents = 0.0;
  158.     double  LyapDimen;
  159.  
  160.     for(i = 0; (SumOfExponents += lyapsum[i]) >= 0.0
  161.             && i < num_lyap; i++);/* sets i */
  162.     if(SumOfExponents < 0) {/* the sum of i exponents is >= 0 *//* notice
  163.                    lyapsum[i] must be negative */
  164.         SumOfExponents = SumOfExponents - lyapsum[i];
  165.         LyapDimen = i + SumOfExponents / (-lyapsum[i]);
  166.         fprintf(output, "Lyapunov Dimension = %lf.  ", LyapDimen);
  167.     }
  168.     else {
  169.         if(num_lyap == vec_dim)
  170.             fprintf(output, "Lyapunov dimension = %d    ", vec_dim);
  171.     }
  172. }                /* see also yintrpt for printout */
  173.  
  174. set_y() {            /* y[j] is set for j = 5,...,20   -- for dif eq
  175.                    223 */
  176.     int    i, coord;
  177.  
  178.     for(i = 0; i < num_lyap; i++)
  179.         for(coord = 0; coord < vec_dim; coord++)
  180.             y[lyapzero + coord + vec_dim * i] = vector[coord][i];
  181. }
  182.  
  183. set_vector() {
  184.     int    i, coord;
  185.  
  186.     for(i = 0; i < num_lyap; i++)
  187.         for(coord = 0; coord < vec_dim; coord++)
  188.             vector[coord][i] = y[lyapzero + coord + vec_dim * i];
  189. }
  190.  
  191. orthonormalize() {        /* orthonormalize  the vectors */
  192.     double  norm;
  193.     int     coord;
  194.     int    i;
  195.  
  196.     product = 1.;
  197.  
  198.     for(i = 0; i < num_lyap; i++) {
  199.         if(i > 0)
  200.             orthogonalize(i);
  201.                 /* make ith vector perpendicular to the earlier
  202.                    vectors */
  203.         norm = 0;
  204.         for(coord = 0; coord < vec_dim; coord++)
  205.             norm += vector[coord][i] * vector[coord][i];
  206.  
  207.         norm = sqrt(norm);
  208.         product *= norm;
  209.         for(coord = 0; coord < vec_dim; coord++)
  210.             vector[coord][i] = vector[coord][i] / norm;
  211.  
  212.         if(dot > 0)
  213.             lyapsum[i] = lyapsum[i] + log(norm);
  214.     }            /* ith vector is finished  */
  215.     if(dot > 0)
  216.         lyaptime = lyaptime + lyapstep;
  217.     return;
  218. }
  219.  
  220.  
  221. /* ORTHOGONALIZE - change vector[][i] into vector[][i] so that
  222.  * vector[][i] is perpendicular to previous
  223.  * ones. This new vector[][i] is unnormalized.
  224. */
  225. orthogonalize(i) 
  226. int    i;
  227. {
  228.     int     k, coord;
  229.     double  dot_p();
  230.  /* dot_prod[k] will be the dot product of the k th vector with the ith already
  231.     normalized vector */
  232.     for(k = 0; k < i; k++) {
  233.         dot_prod[k] = dot_p(i, k);
  234.         for(coord = 0; coord < vec_dim; coord++)
  235.             vector[coord][i] = vector[coord][i] - vector[coord][k] * dot_prod[k];
  236.     }            /*  now we have orthogonal but unnormalized
  237.                    vector # i */
  238. }
  239.  
  240.  
  241. double  dot_p(kk, ii)        /* returns dot product as function value */
  242. int     ii,
  243.         kk;
  244. {
  245.     int     crd;
  246.     double    dotprd = 0.0;
  247.  
  248.     for(crd = 0; crd < vec_dim; crd++)
  249.         dotprd += vector[crd][ii] * vector[crd][kk];
  250.     return(dotprd);
  251. }
  252.  
  253. /* Newton Routines */
  254.  
  255. newtinit() {
  256.     int    coord;
  257.  
  258.     dim = lyapzero + vec_dim * vec_dim;/* necessary for rungekutta */
  259.     for(coord = 0; coord < vec_dim * vec_dim; coord++)/* matrix = 0 */
  260.         y[lyapzero + coord] = 0.;
  261.     for(coord = 0; coord < vec_dim; coord++)/* diagonal elements = 1  */
  262.         y[lyapzero + coord + vec_dim * coord] = 1.;
  263. }
  264.  
  265. #define IF1PRINT   if( printer >=1 ) PRINT
  266.  
  267. double  newton(N, quasi)    /* quasi = 1 means use Quasi-Newton; y1[] is
  268.                    used as the initial state; error is returned
  269.                    */
  270. int     N,
  271.         quasi;            /* N is the number of iterates */
  272. {
  273.  
  274.     int     J;
  275.     double  old[VECDIM],
  276.             new[VECDIM],
  277.             dif[VECDIM],
  278.             miss,
  279.             miss2,
  280.             factor;
  281.     double  ystep[VECDIM];
  282.     double  yy[EQUATIONS];
  283.     int     i,
  284.             n_lyap;
  285.  
  286.     erase_line();
  287.     ResetScrnLineAtTop();
  288.  
  289.     InNewton = 1;        /* indicates we are in newton(); this causes
  290.                    suppression of lyapunov() in iterate_map() 
  291.                 */
  292.  
  293.  
  294.  
  295.     setequal(0, 1, eqn1);    /* set y to equal the initialization vector */
  296.     store(yy, y, zeroth);    /* stores  y in yy -- yy is not currently in
  297.                    use anywhere */
  298.     n_lyap = num_lyap;    /* save num_lyap so it can be reset */
  299.     num_lyap = vec_dim;    /* map iterates look at this number */
  300.  
  301.     newtinit();
  302.  
  303.     for(i = 0; i < vec_dim; i++)
  304.         old[i] = y[zeroth + i];
  305.  
  306.  /* *** now iterate and report resulting error *** */
  307.     for(i = 0; i < N; i++)
  308.         (*iteratee) ();
  309.     for(i = 0; i < vec_dim; i++)
  310.         dif[i] = old[i] - y[zeroth + i];
  311.                 /* for 2 dimensional, the two dimensions are
  312.                    zeroth and zeroth+1 */
  313.     miss = 0;
  314.     for(i = 0; i < vec_dim; i++)
  315.         miss += dif[i] * dif[i];
  316.     miss = sqrt(miss);
  317.     EL;
  318.     IF1PRINT "PERIOD OF ORBIT(set by \"NP\") is %d \n", period);
  319.  
  320.     EL;
  321.     IF1PRINT
  322.         "Before Newton step, a periodic orbit is missed by %le     \n"
  323.         ,miss);
  324.  
  325. /***********************************************************
  326.       lyapzero is probably 2 so the jacobian matrix is
  327.         (  y[2+0]  y[2+2]  )
  328.         (  y[2+1]  y[2+3]  )
  329.     then the inverse of DF-Id is
  330.         (  y[2+3]-1  -y[2+2]    )
  331.         ( -y[2+1]     y[2+0]-1  )     divided by its determinant
  332. ************************************************************/
  333.  
  334. /******** now compute and report recommended newton step(ystep) ************/
  335.     if(vec_dim == 2 && printer > 0)
  336.         newt2 (ystep, N, dif);
  337.     else
  338.         newtMany(ystep, dif);
  339.     EL;
  340.     IF1PRINT
  341.         "Newton step has norm %le \n"
  342.         ,sqrt(ystep[0] * ystep[0] + ystep[1] * ystep[1]));
  343.  
  344.     factor = 1.;
  345.  
  346. /************* Now check to see if we get an improvement in error ************/
  347.     num_lyap = 0;        /* this is so the next iterate does not change
  348.                    the tangent vectors; hence when y[] is
  349.                    displayed using 300, it will display matrix
  350.                    elements */
  351.     dim = lyapzero;        /* rungekutta uses dim */
  352.     setequal(0, 1, zeroth);/* set y to equal the initialization vector;
  353.                    for maps this line can be deleted but for
  354.                    dif eqns it resets time  */
  355.     for(i = 0; i < vec_dim; i++)
  356.         y[zeroth + i] = old[i] + ystep[i] * factor;
  357.     for(i = 0; i < vec_dim; i++)
  358.         new[i] = y[zeroth + i];
  359.  
  360.     setequal(1, 0, zeroth);
  361.     for(i = 0; i < N; i++)
  362.         (*iteratee) ();
  363.  
  364.     for(i = 0; i < vec_dim; i++)
  365.         dif[i] = new[i] - y[zeroth + i];
  366.     miss2 = 0;
  367.     for(i = 0; i < vec_dim; i++)
  368.         miss2 += dif[i] * dif[i];
  369.     miss2 = sqrt(miss2);
  370.  
  371.     EL;
  372.     PRINT
  373.         "AFTER Newton step, a periodic orbit is missed by  %le \n"
  374.         ,miss2);
  375.  
  376.     EL;
  377.     IF1PRINT "\n");
  378.  
  379. /**** now set y and y1 to equal the result of one Newton step ****/
  380.     setequal(0, 1, zeroth);/* reset y to equal the old vector; affects dif
  381.                    eqns that have time as a variable */
  382.  
  383. /**************************  QUASI-NEWTON ***************************/
  384. /**if the error after the Newton step is too big, cut the size ******/
  385.     for(J = 0; (quasi == 1) && (miss2 > miss) && (miss2 >.000001) && (J < 10);
  386.             J++) {
  387.         factor *= 0.5;
  388.         setequal(0, 1, zeroth);/* reset y to equal the old vector */
  389.         for(i = 0; i < vec_dim; i++)
  390.             y[zeroth + i] = old[i] + ystep[i] * factor;
  391.         for(i = 0; i < vec_dim; i++)
  392.             new[i] = y[zeroth + i];
  393.  
  394.         for(i = 0; i < N; i++)
  395.             (*iteratee) ();
  396.  
  397.         for(i = 0; i < vec_dim; i++)
  398.             dif[i] = new[i] - y[zeroth + i];
  399.         miss2 = 0;
  400.         for(i = 0; i < vec_dim; i++)
  401.             miss2 += dif[i] * dif[i];
  402.         miss2 = sqrt(miss2);
  403.  
  404.         EL;
  405.         PRINT
  406.             "Using Newton step TIMES%le, a periodic orbit is missed by %le\n"
  407.             ,factor, miss2);
  408.     }
  409.     setequal(0, 1, zeroth);/* reset y to equal the old vector plus the
  410.                    step times factor  */
  411.     for(i = 0; i < vec_dim; i++)
  412.         y[zeroth + i] = old[i] + ystep[i] * factor;
  413.  
  414.     num_lyap = n_lyap;
  415.     dim = lyapzero + vec_dim * num_lyap;
  416.                 /*  this is needed for "300" for printing out
  417.                    the right number of coordinates of y[] and
  418.                    for the differential equation solver */
  419.     EL;
  420.     IF1PRINT "\n");
  421.     EL;
  422.     IF1PRINT "\n");
  423.     setequal(1, 0, eqn1);    /* set the initialization vector equal to y */
  424.  
  425.     InNewton = 0;        /* indicates newton() is being left */
  426.     return(miss2);
  427. }
  428.  
  429. newtMany(ystep, dif)        /*  mat = df - id,  mat*ystep = dif = old -
  430.                    (y+zeroth) */
  431. double *ystep,
  432.        *dif;
  433. {
  434.     int     i,
  435.             j;
  436.     double  mat[MATDIM];
  437.  
  438.     for(i = 0; i < vec_dim; i++)/* for gauss() */
  439.         for(j = 0; j < vec_dim; j++)
  440.             mat[i + j * vec_dim] = y[lyapzero + i * vec_dim + j];
  441.     for(i = 0; i < vec_dim; i++)/* subtract identity */
  442.         mat[i + i * vec_dim] -= 1.0;
  443.     gauss(mat, ystep, dif, vec_dim);
  444.                 /* solve mat x = b where vec_dim is the
  445.                    dimension of the vectors and if for vec_dim
  446.                    = 2, mat = ( mat[0] mat[1] ) ( mat[2] mat[3]
  447.                    ) WARNING: mat is altered though b is not */
  448.     for(i = 0; i < vec_dim; i++)/* for gauss() */
  449.         EL;
  450.     IF1PRINT
  451.         "after;[%d]:ystep=%lf dif=%lf\n", i, ystep[i], dif[i]);
  452.  
  453.     return;
  454. }
  455.  
  456.  
  457. newt2 (ystep, N, dif)        /* computes ystep[] in two dimensions and
  458.                    prints out eigenvalues,etc */
  459. double *ystep,
  460.        *dif;
  461. int     N;
  462. {
  463.     double  A,
  464.             B,
  465.             C,
  466.             D,
  467.             det;
  468.     double  a,
  469.             b,
  470.             c,
  471.             d;
  472.  
  473.     a = y[lyapzero];
  474.     b = y[lyapzero + 2];
  475.     c = y[lyapzero + 1];
  476.     d = y[lyapzero + 3];
  477.     EL;
  478.     IF1PRINT
  479.         "The Jacobian Matrix for  %d iterate(s) is\n", N);
  480.     EL;
  481.     IF1PRINT
  482.         "    ( %lf, %lf )\n", a, b);
  483.     EL;
  484.     IF1PRINT
  485.         "    ( %lf, %lf )     Jacobian det = %lf\n", c, d, a * d - b * c);
  486.  
  487.     radical = (a + d) * (a + d) - 4 * (a * d - b * c);
  488.     EL;
  489.     if(radical > 0) {    /* distinct real eigenvalues */
  490.         EigenValue1 = (a + d - sqrt(radical)) / 2;
  491.         EigenValue2 = (a + d + sqrt(radical)) / 2;
  492.         EL;
  493.         IF1PRINT
  494.             "The eigenvalues are     %lf   and   %lf \n"
  495.             ,EigenValue1, EigenValue2);
  496.         if(b != 0) {
  497.             Ycomponent1 = -(a - EigenValue1) / b;
  498.             Ycomponent2 = -(a - EigenValue2) / b;
  499.             EL;
  500.             IF1PRINT
  501.                 "    with eigenvectors:  (1, %lf)    (1, %lf)    \n",
  502.                 Ycomponent1, Ycomponent2);
  503.         }
  504.     }
  505.     if(b == 0 && c != 0) {
  506.         EL;
  507.         IF1PRINT
  508.             "with eigenvectors:      (%lf, 1)    (%lf, 1)\n",
  509.             -(d - EigenValue1) / c, -(d - EigenValue2) / c);
  510.     }
  511.     if(radical < 0) {    /* complex eigenvalues */
  512.         EL;
  513.         IF1PRINT
  514.             "The eigenvalues are complex and = %lf +- i%lf \n"
  515.             ,(a + d) / 2, sqrt(-radical) / 2);
  516.     }
  517.  
  518.  
  519.     A = y[lyapzero + 3] - 1;
  520.     B = -y[lyapzero + 2];
  521.     C = -y[lyapzero + 1];
  522.     D = y[lyapzero] - 1;
  523.     det = A * D - B * C;
  524.  
  525.     if(det <.00000000001 && det > -.00000000001) {
  526.         EL;
  527.         IF1PRINT
  528.             "+1 is an eigenvalue; cannot take Newton step\n");
  529.         InNewton = 0;
  530.         return;
  531.     }
  532.     EL;
  533.     IF1PRINT
  534.         "The Newton Matrix, i.e. the inverse of Jacobian-Id, is \n");
  535.  
  536.     A = A / det;
  537.     B = B / det;
  538.     C = C / det;
  539.     D = D / det;
  540.  
  541.     EL;
  542.     IF1PRINT
  543.         "    ( %lf, %lf )\n", A, B);
  544.     EL;
  545.     IF1PRINT
  546.         "    ( %lf, %lf )   with determinant %lf\n", C, D, 1 / det);
  547.  
  548.  
  549. /*    (A,B,C,D) = inverse of [Jacobian - Identity]    */
  550. /*    y = old + (A,B,C,D)*(old - y)            */
  551.  
  552.     ystep[0] = A * dif[0] + B * dif[1];
  553.     ystep[1] = C * dif[0] + D * dif[1];
  554.  
  555.     return;
  556. }
  557.  
  558. /* %* %* %* %* %* %* %* %* %* %* %* %* %* %* %* %* %* %* %* %* */
  559. /* Routine gauss() is to use gaussian method to solve system   */
  560. /* of linear equations. It contains three subroutines which are*/
  561. /* pivot(),tri() and solve(). Pivot() is used to pivot each    */
  562. /* coming column, tri() is to be used to transform the matrix  */
  563. /* into upper triangular matrix, and solve() is to use chasing */
  564. /* method to solve the linear system with an upper triangular  */
  565. /* matrix on the left hand side.                           */
  566. /*                                        */
  567. /* Notice: We input all the matries as vectors, i.e. an N by N */
  568. /* matrix is to be considered as a vector with dimension N*N,  */
  569. /* and we order them row by row. Therefore the element a[i][j] */
  570. /* of an ordinary matrix is denoted as a[i*vec_dim+j] here.    */
  571. /*                                                             */
  572. /*                                                             */
  573. /*  To use this subroutine you should first define             */
  574. /* VECDIM & MATDIM, where VECDIM = dimension of the vector     */
  575. /* (the number of variables) and MATDIM = the dimension of the */
  576. /* matrix, i.e. MATDIM=VECDIM*VECDIM. The format of define is  */
  577. /* #define  VECDIM  3 . Then you should give                */
  578. /* values to the a[i]'s. Now you simply call routine gauss()   */
  579. /* to solve your linear equation. Here is an example:          */
  580. /*                                                             */
  581. /*     #include <stdio.h>                                      */
  582. /*                                                             */
  583. /*     #define  VECDIM  4                                      */
  584. /*     #define  MATDIM  16                                      */
  585. /*                                                             */
  586. /*     solve  ax = b                                           */
  587. /*     main()                                                  */
  588. /*     {                                                       */
  589. /*           double a[MATDIM],b[VECDIM],x[VECDIM];             */
  590. /*                                                             */
  591. /*           a[0]=2;                                           */
  592. /*           a[1]=0;                                           */
  593. /*           a[2]=1;                                           */
  594. /*           a[3]=1;                                           */
  595. /*           b[0]=0;                                           */
  596. /*           b[1]=1;                                           */
  597. /*                                                             */
  598. /*           gauss(a,x,b);                                     */
  599. /*           printf("%lf,%lf\n",x[0],x[1]);                    */
  600. /*      }                                                      */
  601. /* %* %* %* %* %* %* %* %* %* %* %* %* %* %* %* %* %* %* %* %* */
  602.  
  603.  
  604. gauss(a, x, b, vec_dim)    /* solve ax = b where vec_dim is the dimension
  605.                    of x and b and where for vec_dim = 2, a = (
  606.                    a[0] a[1] ) ( a[2] a[3] ) WARNING: a is
  607.                    altered though b is not */
  608. double *a,
  609.        *x,
  610.        *b;
  611. int     vec_dim;
  612. {
  613.     int     j;
  614.  
  615.     for(j = 0; j < vec_dim; j++) {
  616.         x[j] = b[j];    /* now b will not be changed */
  617.     }
  618.  
  619.     for(j = 0; j < vec_dim; j++) {
  620.         pivot(a, x, j, vec_dim);
  621.         tri(a, x, j, vec_dim);
  622.     }
  623.     solve(a, x, vec_dim);
  624.     return;
  625. }
  626.  
  627. /* %* %* %* %* %* %* %* %* %* %* %* %* %* %* %* %* %* %* %* %* */
  628. /* subroutine pivot() pivots the j-th column of the current matrix*/
  629. /* a[i][j] (with i>=j) so that in tri() we can avoid to devide */
  630. /* something by zero. Furthermore, because of pivoting, we can */
  631. /* (in some extent) minimize the accumutive error due to round */
  632. /* off by computer.                                            */
  633. /* %* %* %* %* %* %* %* %* %* %* %* %* %* %* %* %* %* %* %* %* */
  634.  
  635.  
  636. pivot(a, z, j, vec_dim)
  637. double *a,
  638.        *z;
  639. int     j,
  640.         vec_dim;
  641. {
  642.     int     i,
  643.             l,
  644.             k,
  645.             jj,
  646.             kk;
  647.     double  s[MATDIM],
  648.             x[VECDIM],
  649.             b[VECDIM],
  650.             d,
  651.             aa;
  652.  
  653.  /* s and x are just passing variables. When we switch two    */
  654.  /* rows, we put the first row into the passing variables and */
  655.  /* the put the second row into the first one, and finally we */
  656.  /* put the the passing variables' value into the second row. */
  657.  
  658.     for(i = j; i < vec_dim; i++) {
  659.         aa = a[i * vec_dim + j];
  660.         (aa < 0) ? (b[i] = -aa) : (b[i] = aa);
  661.     }
  662.  
  663.  /* b is the absolute value of the coresponding a[i][j] so to */
  664.  /* make the following comparison possible.(To find the one   */
  665.  /* with the largest absolute value in the j-th column).      */
  666.  
  667.     k = j;
  668.     d = b[j];
  669.     for(i = j + 1; i < vec_dim; i++)
  670.         if(b[i] > d) {
  671.             d = b[i];
  672.             k = i;
  673.         }
  674.     if(d == 0) {
  675.         printf("pivot(): The system is unsolvable because the \
  676. matrix is singular!");
  677.         exit(1);
  678.     }
  679.  /* the following is to switch the entries of the rows due to */
  680.  /* partial pivoting.                         */
  681.     for(l = j; l < vec_dim; l++) {
  682.         jj = j * vec_dim + l;
  683.         kk = k * vec_dim + l;
  684.         s[jj] = a[jj];
  685.         a[jj] = a[kk];
  686.         a[kk] = s[jj];
  687.     }
  688.     x[j] = z[j];
  689.     z[j] = z[k];
  690.     z[k] = x[j];
  691.     return;
  692. }
  693.  
  694. /* %* %* %* %* %* %* %* %* %* %* %* %* %* %* %* %* %* %* %* %* */
  695. /* subroutine tri() is to transform the matrix into an upper   */
  696. /* triangular one so that it is easy to solve. What we really  */
  697. /* do in this subroutine is to transform the j-th column of the*/
  698. /* current matrix into the form such that all elements a[i][j] */
  699. /* with i>j are zero.                                          */
  700. /* %* %* %* %* %* %* %* %* %* %* %* %* %* %* %* %* %* %* %* %* */
  701.  
  702.  
  703. tri(a, x, j, vec_dim)
  704. double *a,
  705.        *x;
  706. int     j,
  707.         vec_dim;
  708. {
  709.     int     i,
  710.             k;
  711.  
  712.     for(i = j + 1; i < vec_dim; i++) {
  713.         for(k = j + 1; k < vec_dim; k++)
  714.             a[i * vec_dim + k] = a[i * vec_dim + k]
  715.                 - a[j * vec_dim + k] * (a[i * vec_dim + j] / a[j * vec_dim + j]);
  716.         x[i] = x[i] - x[j] * (a[i * vec_dim + j] / a[j * vec_dim + j]);
  717.     }
  718.     return;
  719. }
  720.  
  721.  
  722. /* %* %* %* %* %* %* %* %* %* %* %* %* %* %* %* %* %* %* %* %*   */
  723. /* Subroutine solve() is to use chasing method to solve a system */
  724. /* with an upper triangular matrix on the left hand side. We     */
  725. /* solve the n-th unknown first, and put it back in x[n], and    */
  726. /* solve n-1th unknown and so on.                                */
  727. /* %* %* %* %* %* %* %* %* %* %* %* %* %* %* %* %* %* %* %* %*   */
  728.  
  729.  
  730. solve(a, x, vec_dim)
  731. double *a,
  732.        *x;
  733. {
  734.     int     i,
  735.             k;
  736.  
  737.     for(i = vec_dim - 1; i > -1; i--) {
  738.         for(k = vec_dim - 1; k > i; k--)
  739.             x[i] -= a[i * vec_dim + k] * x[k];
  740.         x[i] /= a[i * vec_dim + i];
  741.     }
  742.     return;
  743. }
  744.